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Abstract We developed a least squares fitter used for extracting expected physics parameters from the 
correlated experimental data in high energy physics. This fitter considers the correlations among the observables 
and handles the nonlinearity using linearization during the x^ minimization. This method can naturally be 
extended to the analysis with external inputs. By incorporating with Langrange multipliers, the fitter includes 
constraints among the measured observables and the parameters of interest. We applied this fitter to the study 
of the D" — -D" mixing parameters as the test-bed based on MC simulation. The test results show that the 
fitter gives unbiased estimators with correct uncertainties and the approach is credible. 
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1 Introduction 



It frequently happens that one wants to deter- 
mine the unknown parameters from a set of correlated 
experimental measurements. Least squares fit [1] is 
an effective and standard approach for this purpose. 
The most general situation is the estimation problem 
involving the observables and unknown parameters, 
which are connected through a set of linear and non- 
linear constraints. It is well known that if the con- 
straints are linear equations, least squares fit gives 
unbiased results with correct uncertainties. For non- 
linear constraints, minimization becomes more com- 
plex and linearization are often introduced so that 
it can be solved by linear solutions. However, those 
results from linearization can be slightly biased and 
variances are not minimum in general. Thus, good 
approximation in the linearization is required. 

For data analysis in high energy physics experi- 
ments, the observables are mostly event numbers and 
their relations with the parameters of interest are 
nonlinear in most cases. Furthermore, global fit is 
an important method to better constrain the param- 
eters by combining the experimental measurements 
and the external inputs. In this paper, we develop 



an approach based on least squares fit and Langrange 
multiplier method for these cases. The statistical and 
systematic uncertainties of the observables and their 
dependencies on the fit parameters [2| are considered 
in constructing the characteristic x^ ^'^d the mini- 
mization procedure. 

2 Formalism 

Throughout this paper, the lowercase bold letter 
refers to vector quantity, the uppercase letter repre- 
sents matrix quantity, the symbol V stands for co- 
variance matrix. 

2.1 Construction of x^ 

In least squares fit with constraints, the unknown 
parameters m can be obtained by minimizing x^- Re- 
ferring to Ref [3|, |j], we construct the x^ in an ex- 
tended form: 



X' = (y-»7rV^'(y-r7) + 2A„^g(r7,m) 
+ 2A/h(r,), 



(1) 



where y is the vector of experimental observations, 
and T] is the expected value of y. Generally, jy is a 
function of m and their relationship can be expressed 
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as g(j7, m)=0. h{ri) is the vector of constrain func- 
tions of T]. Xd and A^ are the vectors of Langrange 
multipliers. Minimizing the x^ leads to find the opti- 
mized value of m. Typically, Vy is determined from 
experimental measurements and is taken as a con- 
stant in x^ fit- However, there are cases that Vy 
depends on m. With different input m, the weight of 
each measurement should be altered. Otherwise, the 
result may be biased. In our case, we consider Vy 
as Vy(m), and it will be updated in each iteration of 
the fit. 

Let's discuss the usual cases of measurements in 
high energy physics experiments, where direct observ- 
ables are the numbers of signal events n. Each item in 
n corresponds to the number of event candidates of a 
physics process. With extraction of the backgrounds, 
their expected values are functions of m, which in 
most cases are branching fractions. Usually, the sig- 
nal events n may receive crossfeed contributions from 
other signal processes and contaminations from peak- 
ing backgrounds which are not belonging to the pro- 
cesses of interest. We use b to describe the number of 
these peaking backgrounds. The efficiencies-corrected 
yields, denoted by c, can be expressed as: 

c = E"'s = E"'(n-Fb), (2) 

where, E is the signal efficiencies matrix, to describe 
detection efficiencies and crossfeed probabilities, F is 
background efficiencies matrix, to describe contami- 
nation rates from background to each signal process. 
Assuming that there are external measurements 
t that can be incorporated to constrain parameters 
of interest further, the x^ can be built with all the 
measurements c and t included in y: 



ti 



(3) 



In the case of that g(j7, m) is nonlinear, Taylor 
expansion to the first order can be given as: 

Qg 9g 

g(?7,m)«g(rjo,mo)+^-— (m-mo)+T^(r7-r7(,). (4) 
om or] 

Here we assume that the deviation from point (m, r/) 
to (mo, J7o) should be small. The similar linearization 
is also applied on h(j7). 

2.2 Input variance 

To obtain unbiased fit results, proper handling of 
variance matrixes is required. According to Eq.(2), 



the uncertainties of n, b, E, F should be propagated 
to c as: 



V =(^VY ^+(^rv,— 




(5) 

where, V„, V^, Vg, Vp are the uncertainties of n, 
b, E, F respectively. Generally, the variances of E 
and F depend on uncertainties of estimating track- 
ing efficiency, particle identification(PID) and so on. 
For the different processes, the uncertainties of the 
observables are correlated. Therefore, Ve and Vp 
have nonzero off-diagonal elements. Also E and F 
share many common correlated uncertainties. These 
common uncertainties are denoted by Cef- More dis- 
cussions about Vc can be found in Ref. [2|. 

In general cases, external measurements are not 
related to the internal measurements. Therefore, Vy 
is simplified as: 



V. 




(6) 



Vt is the variance matrix of t. In the case of correla- 
tion exist between c and t, the off-diagonal elements 
should be nonzero. 

2.3 Minimizing x^ 

There are many approaches in x^ minimization. 
We adopt the iterative procedure. That is, the esti- 
mated values in step fc, m'°, are used as seeds for cal- 
culating the estimators m*+'^ in the step k + 1. The 
equation is formulated as [3|, |j] : 



m'=+^ =m^ - [GlS-^GirrGlS-^ 
[zi-(G^)'=Vy(H^)^-S-^Z2], 
where 

Si = (G^)'=VyG:;, 

s,^(H^)'=VyH:; 

S3^(G^)'=VyH^S-(H^)'=Vy(G;), 
S4 = Si — S3, 

Zi^g^ + G^(y-r,'=), 
Z2EEh'=+H^(y-,7^-)- 



(7) 



(8) 

(9) 

(10) 
(11) 
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The fit procedure is to reiterate Eq.(7) until the x^ 
converges. Then we obtain the variance matirx as 



-S.Y^Si, 



where 



S5=[GmS4 Gjjj] 0^84 G^ 
[I-V^H^S-^H^], 



(12) 



(13) 



where I indicates the unit matrix. More details about 
deducing Eq.(7-13) are put in the Appendix A. In our 
specific case, c is dependent on m(through b). Note 
that dc/ dva can be ignored in x^ minimization, be- 
cause the elements of F are very small in general. 
9Vc/9m is not considered in deriving Eq.(7). This 
special treatment avoids the potential bias [2] , which 
is introduced by this item. However, in each iteration 
all the input variables that depend on m are recalcu- 
lated, including Vc and c. 

3 Monte Carlo study 

The fitter is developed based on ROOT [5| frame- 
work. We test it in the measurement of D° — Z)" 
mixing parameters by toy Monte Carlo (MC) simu- 
lation under the environment of the BESIII experi- 
ment [6|, where D-pair is produced through e+e~ -^ 
^^(3770) — )► DD, and they are in a quantum-correlated 
C-odd system |7|, [Sf . The measurement of their de- 
cay rates provide unique opportunity for measuring 
£)0_£)0 ]-|2ixing parameters [9|4ll|. We use ten signal 
processes as listed in Table [T] [l^ . 

Table 1. Signal processes involved in the test. 
/""'■ are the correlated (C-odd) effective D^D'^ 
branching ratios, to the leading order in xd, 
yo and Rws, divided by the branching ratios 
Bi of a isolated D for modes i and BiBj for 
modes {i, j}. 



D decay mode 




£cor 


X~7r+ 




l + Rws 


K+K- 




2 


KsttO 




2 


K~n+,K+n'' 


(1 + Rwsf- 


- 4r cos Sxvrir cos Sktv + VD ) 


K-7T+,K+K- 


1 + Rws + 2»- cos <5x,r + Vd 


K-n+,Ks7rO 


1 + Rws - Sr cos Sji^ - yo 


K-TT+,K+e-pe 


^-ryo 


cosSxTT-rxosinSKn 


K+K-,KsnO 




4 


K+K-,KeUe 




2(1 + S/I3) 


Ksn°,KeUe 




2(1-S/d) 



^D, Vd- Nod indicates the total number of produced 
D^D^ pairs; B indicates the branching ratios; —Skw 
is the relative phase between the doubly Cabibbo- 
suppressed D° — >■ isT+Tr" amplitude and the corre- 
sponding Cabibbo-favored D" -^ K'^tt^ amplitude: 
<K+TT~\D''> / <K+n-\D''>=re-''^'^^] Xd, Vd are 
parameters describes charm mixing, for the details 
of these definitions, we refer to Ref [1^. We input 
NuD = 5.0 X 10'', which roughly corresponds to those 
yields in S.O/fe-^ data of e+e' ^DD at the V(3770) 
resonance. The values of other input parameters are 
taken as the world-average values [13|] with Gaussian 
smearing. The width of Gaussian is taken as the er- 
ror of the corresponding parameter. Detection effi- 
ciencies for these processes are determined from MC 
sample of simulating the BESIII detector. We as- 
sume 0.5% peaking backgrounds (from p7r processes) 
for the modes involved with D -^ KsTr°. We apply 
correlated systematic uncertainties of 1% for track- 
ing efficiencies, 2% for tt" finding and 4% for Ks find- 
ing. All the event yields are fluctuated according to 
Poisson statistics. 

In the fit to the MC sample, we take inputs of 
data from other experimental measurements, which 
can provide more constraints on parameters of in- 
terest. There are seven external inputs in the test: 
Rwsi f^, Sktv, Xo, 2/d, x'^ and y' and their uncertain- 
ties are assumed to be uncorrelated. Relationships 
among these variables are listed as follows: 



Rws =r'^ + ryo cos{5Kn) - rxo sin((5K^) 
+ 0.5x(4+yJ), 

X' ~Xo cos 5 Kit + J/D Sin Sk-k , 

y' ^yDCOs5KT,-XDSin5K^. 



(14) 



The fit is expected to reproduce nine parameters: 
Nod, BiKn), B{KK), S(i^s^°), B{Keu), r, <5^,, 



We do ten thousands times of sampling and per- 
form the least squares fit for each sample. The pull 
distributions for nine fit parameters are shown in 
Fig. [1] All the pull distributions agree well with 
the normal distribution and the confidence level is 
fiat. This indicates that the fitter provides unbiased 
estimations of the parameters of interest and good 
convergence. Slight asymmetries in pull distributions 
may present, due to the nonlinearity. Table [5] lists the 
correlation coefficients among the fit parameters. As 
we expect, branching fractions tend to be positively 
correlated with each other and negatively correlated 
with Ndd- 
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Fig. 1. Pull distributions of Ndd{^), 
B{K-K){h), B{KK){c), S(is:s7r°)(d), 

B{Keu){e), r(f), 5K„{g), a;i3(h), yoii), 
overlaid with normal distributions. The 
confidence level distribution (j) is overlaid 
with a line with zero slope. 

We also estimate the sensitivity of measuring 
Hd and Sktv under the current statistics. Consider- 
ing more available modes, in this estimation, events 



yields for CP eigenstates and semi-leptonic processes 
are scaled by a factor of 2 roughly. We input world- 
average Sk^ = 22.1+?i';i(°) and yo = 0.75±0.12(%) ^ 
for the fit test. One-dimensional confidence curves of 
the fit of yo and Sktt are shown in Fig. [2j The curves 
are obtained by repeating the fits at fixed value of y^ 
or Sktt and recording the change form the minimum 
Xmin- The uncertainties of output Sktv and yo are 
determined to be 8.3° and 0.10% respectively, which 
are determined by Ax^ = 1. The results show that 
uncertainties on y^ and (5x^ are both improved by 
about 15%. 




Fig. 2. The function Ax^=X^-Xmin for yo(a) 
and 5Kir(b). The dashed line denotes the 
points where Ax^=l. 



Table 2. Correlation coefficients, including systematic uncertainties, for the parameters determined by the fit 
with MC samples. 



Ndd 


B{Ktt) 


B{KK) 


B{Ks-K°) 


B{Keu) 
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xo 
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-0.01 


-0.01 


B{K-k) 


1 


0.96 


0.15 


0.65 


0.01 


-0.02 


0.00 


0.01 


B{KK) 
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0.15 
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0.00 


0.01 


B{KsnO) 
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0.01 
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0.00 


0.01 


B{Keu) 
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-0.00 


0.00 


-0.00 


0.01 
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-0.28 


Sk. 
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0.09 


xn 














1 


-0.09 


Vd 
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4 Summary 

We developed a least squares fitter, which extracts 
the expected parameters by combining the experi- 
mental measurements and the external inputs. Lan- 
grange multiplier method is adopted accounting for 
constraints among the observables and the expected 
parameters. In the fitter, the observables and the in- 



put covariance matrix are supposed to be dependent 
with the expected parameters and they need to be 
renewed in each iteration step during minimization 
procedure. With correct input of the error matrix of 
the observables, the fitter gives unbiased estimations 
with correct uncertainties of the expected parame- 
ters. The test on toy MC validates the credibility of 
the fitter. 
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Appendices A 

Formulas for iterative process 

Following the similar procedure presented in Ref [3,1^, 
one can obtain Eq.(7-13). By assuming the deviation of 
~X^ to 7j, m, Aa, A(3 equal to zero, we obtain 

-2V;'(y-77) + 2G^A„ + 2H^A0 = O, (Al) 

2G^A„=0, (A2) 

2g(7,,m) = 0, (A3) 

2h(ry) = 0. (A4) 

m*"*"^, T/*"*"^, AJ^+^, and A^+^ are used as inputs to the 
next iteration. Eq.(Al) and (A2) can be re-expressed as: 



V-(r7'=+^- 



-y) + G:;A 



fc+i 



ttA: ^ fc+l 



= 0, 



(A5) 

C = 0. (A6) 

With Taylor expansion, Eq.(A3) and (A4) become: 



g^ + (G^)''-(77'=+^-r,'=) + (G^)'=(m'=+i-m'=)=0, (A7) 

h'= + (H^)'=(ry'=+'-ry'=)=0. 
From Eq.(A5), we have 



V 



(A8) 
(A9) 



y VyGjjAa VyH^^A^ 

With input of Eq.(A9), Eq.(A7) and Eq.(A8) are re- 
written as: 

(AlO) 



zi — SiAq, — (Gj,) VyHj,A^ 



+ (G4,)'=(m'=+^-m'=) = 



Z2 - (H^ )^ VyG:,Ar ' - S2 A;+^ = 0, 



Then 



k ^ fe + l\ 



A;+^ = S2-^(z2-(H^)*^VyG^A; 

A;^+^ in Eq.(AlO) is substituted as 

zi - S4A^+' - (G^)'=Vy(H^)'=S2-iz2 
+ (G^)'=(m'=+^-m''-)=0. 



Then AJL+^ becomes 



(All) 



(A12) 



(A13) 



,fe+i , 



:Sr^[zi-(G^)*=Vy(H^)'=S2-iz2 



(G^) (m ^ -m )]. 



(A14) 



Combing Eq.(A14) and Eq.(A6), we would drive out 

yfc4 



m'=+^ in Eq.(7). The estimators ■q''^^ and A^+^ are ob- 



tained from Eq.(A9) and Eq.(A12). 
The variance matrixes Vm and 
lated variances can be obtained from Eq.(7) and Eq.(A9): 

(A15) 






V -f^rv A 



covir),m)^i^fVy{^) 



^dy 



dy 



(A16) 

(A17) 



References 

1 G. H. Golub, C. F. Van Loan. SIAM Journal on Numerical 
Analysis, 1980, 17: 883-893 

2 W. M. Sun. Nucl. lustrum. Meth. A, 2006, 556(1): 325-330. 

3 Yongsheng Zhu. Probability and Statistics in Experimen- 
tal Physics. 2nd edition. Beijing: Academic Press, 2006. 
Section 9.8(in Chinese) 

4 A. G. Frodesen, O. Skjeggestad, H. Tofte. Probability and 
Statistics in Particle Physics. Universitctsforlaget, Bergen- 
Oslo-Tromso, 1979. Section 10.8 



5 |http://root.cern.ch/download/doc/14Linear Algebra.pdf I 

6 D. M. Asner. International Journal of Modern Physics A, 
2009, 24, 1 

7 Haibo Li. Nuclear Physics B, 2006, 162: 312-346 

8 Jonathan L. Rosner. Phys. Rev. D, 1977, 15, 1254 

9 Bin Huang. Chinese Physics C, 2008, 32(12): 945-951 

10 D. M. Asner, W. M. Sun. Phy Rev. D, 2006, 73, 034024 

11 Z. Z. Xing. Phys. Rev. D, 1997, 55, 196 

12 D. M. Asner et al. Phys. Rev. D, 2012, 86, 112001 

13 J. Bcringcr et al. Phys. Rev. D, 2012, 86, 010001 



